# coding: utf-8
import numpy
import numpy as np
import matplotlib.pyplot as plt
from pandas import read_csv
import math
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from keras.layers import RNN, SimpleRNN, SimpleRNNCell
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
from keras.utils.vis_utils import plot_model


# 创建数据集pi
def create_dataset(dataset):
    train_size = int(len(dataset) * 0.75)
    train, test = dataset[0:train_size], dataset[train_size:len(dataset)]

    return numpy.array(train), numpy.array(test)


def mape(y_true, y_pred):
    return numpy.mean(numpy.abs((y_pred - y_true) / y_true)) * 100


def mean_absolute_percentage_error(y_true, y_pred):
    y_true, y_pred = np.array(y_true), np.array(y_pred)
    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100


def get_r2_numpy(x, y):
    slope, intercept = np.polyfit(x, y, 1)
    r_squared = 1 - (sum((y - (slope * x + intercept)) ** 2) / ((len(y) - 1) * np.var(y, ddof=1)))
    return r_squared


def score_R2(ans, y_test):
    pmean = np.mean(y_test)
    omean = np.mean(ans)
    SSR = 0.0
    varp = 0.0
    varo = 0.0
    for i in range(0, len(y_test)):
        diffXXbar = y_test[i] - pmean
        difYYbar = ans[i] - omean
        SSR += (diffXXbar * difYYbar)
        varo += diffXXbar ** 2
        varp += difYYbar ** 2
    SST = math.sqrt(varo * varp)
    return (SSR / SST) ** 2

if __name__ == '__main__':
    # 加载数据
    dataframe = read_csv('GF-Guy_HH_TOTAL.csv')
    len_dat = dataframe.shape[0]
    # dataframe = dataframe[three_one_len_dat:three_one_len_dat * 2]
    # dataframe = dataframe[three_one_len_dat * 2:len_dat]

    datasetX = dataframe.loc[:, [x for x in dataframe.columns.tolist() if x != 'NEE']].as_matrix()
    datasetY = dataframe.loc[:, [x for x in dataframe.columns.tolist() if x == 'NEE']].as_matrix()

    scaler_X = MinMaxScaler(feature_range=(0, 1))
    scaler_Y = MinMaxScaler(feature_range=(0, 1))
    datasetX = scaler_X.fit_transform(datasetX)
    datasetY = scaler_Y.fit_transform(datasetY)

    # 将整型变为float
    datasetX = datasetX.astype('float32')
    datasetY = datasetY.astype('float32')

    trainX, testX = create_dataset(datasetX)
    trainY, testY = create_dataset(datasetY)

    # 调整输入数据的格式
    # [samples, timesteps, features]
    # trainX = trainX.reshape((trainX.shape[0], trainX.shape[1], 1))
    # testX = testX.reshape((testX.shape[0], testX.shape[1], 1))
    # 创建LSTM神经网络模型
    model = Sequential()
    model.add(Dense(20, input_dim=trainX.shape[1], activation='relu'))
    model.add(Dense(1))
    # model.add(Dense(1, input_dim=trainX.shape[1], activation='relu'))
    # model.add(Dense(1, activation='relu'))
    # model.add(Dense(1, activation='sigmoid'))
    model.compile(loss='mean_squared_error', optimizer='adam', metrics=['accuracy'])
    model.fit(trainX, trainY, epochs=500, batch_size=400, verbose=2)

    # 预测
    trainPredict = model.predict(trainX)
    testPredict = model.predict(testX)

    y_train_true = scaler_Y.inverse_transform(trainY)
    y_pred_true = scaler_Y.inverse_transform(testY)
    y_train = scaler_Y.inverse_transform(trainPredict)
    y_pred = scaler_Y.inverse_transform(testPredict)

    print("train rmse: ", numpy.sqrt(mean_squared_error(y_train, y_train_true)))
    print("train mae: ", mean_absolute_error(y_train, y_train_true))
    print("train mape: ", mean_absolute_percentage_error(y_train_true, y_train))
    print("train r2: ", score_R2(y_train, y_train_true))

    print("predict rmse: ", numpy.sqrt(mean_squared_error(y_pred, y_pred_true)))
    print("predict mae: ", mean_absolute_error(y_pred, y_pred_true))
    print("predict mape: ", mean_absolute_percentage_error(y_pred_true, y_pred))
    print("predict r2: ", score_R2(y_pred, y_pred_true))

    true_list = numpy.array(numpy.concatenate((y_train_true, y_pred_true)))
    pred_list = numpy.array(numpy.concatenate((y_train, y_pred)))

    print("true_list: " + str(list(true_list.squeeze())))
    print("pred_list: " + str(list(pred_list.squeeze())))

    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    ax.plot(range(len(true_list)), true_list, c="black")
    ax.plot(range(len(pred_list)), pred_list, c="red")
    x_ticks = ax.set_xticks(range(len(true_list) / 24, len(true_list), len(true_list) / 12))
    x_labels = ax.set_xticklabels(range(1, 13, 1), fontsize="small")
    ax.set_xlabel("Month")
    ax.set_ylabel("NEE(gCO2 m-2 d-1)")
    ax.legend()
    plt.show()
